---
title: "Replication_testingH1"
author: "Myung Jung Kim"
date: "2026-01-30"
output: html_document
---

```{r R packages setup, include=FALSE}
library(readr)
library(tidyverse)
library(marginaleffects)
library(sandwich)
library(graphics)
library(stats)
library(ggplot2)
library(ggeffects)
library(sandwich)
library(lmtest)
library(stargazer)
library(texreg)
```

```{r H1: Split Models}
df2 <- read_csv("df1_dyadlevel_amnesty.csv")

################################################################
################### A: all data (assume all commit SV)  ##########
################################################################

dfpre98war <- df2 %>%filter(end_pre1998==1)
dfpost98war <- df2 %>%filter(end_pre1998==0)

## ################## Pre98: AX TBR_NSAdummy 
xf1A_98_w_controls_pre89df_AddrebOSV <- glm(hram ~ TBR_NSAdummy
                       + numdyads +rebcap+ blood
                       + Home_v2svstterr+v2x_rule_Home
                       +intensitylevel +dyad_yearatwarln  
                       +PKO_dummy_Home +civ_deaths_dummy 
                      , data=dfpre98war, family = binomial)

# Compute robust standard errors and store the coeftest results
xf1A_98_w_controls_pre89df_AddrebOSV_cluster <- coeftest(
  xf1A_98_w_controls_pre89df_AddrebOSV, 
  vcov = sandwich::vcovCL(xf1A_98_w_controls_pre89df_AddrebOSV, cluster = dfpre98war$side_a_id, type = "HC0")
)
######### Post98: AX TBR_NSAdummy   #####################
xf1A_98_w_controls_post89df_AddrebOSV <- glm(hram ~ TBR_NSAdummy
                       + numdyads +rebcap+ blood
                       + Home_v2svstterr+v2x_rule_Home
                       +intensitylevel +dyad_yearatwarln
                       + +PKO_dummy_Home 
                      +civ_deaths_dummy 
                      , data=dfpost98war, family = binomial)

# Compute robust standard errors and store the coeftest results
xf1A_98_w_controls_post89df_AddrebOSV_cluster <- coeftest(
  xf1A_98_w_controls_post89df_AddrebOSV, 
  vcov = sandwich::vcovCL(xf1A_98_w_controls_post89df_AddrebOSV, cluster = dfpost98war$side_a_id, type = "HC0")
)

```

```{r H1: Interaction Models}
############################################################
######## A: all data   ##########
##############################################

########## [WO control] AX TBR_NSAdummy  X 1998 ##############
xf1A_98_WO_controls <- glm(hram ~ war_ongoing98*TBR_NSAdummy
                   , data=df2, family = binomial)

# Compute robust standard errors and store the coeftest results
xf1A_98_WO_controls_cluster <- coeftest(
  xf1A_98_WO_controls, 
  vcov = sandwich::vcovCL(xf1A_98_WO_controls, 
                          cluster = df2$side_a_id, type = "HC0")
)


## ########### AX TBR_NSAdummy  X 1998 #########################
xf1A_98_w_controls_exceptPKO_AddrebOSV <- glm(hram ~ war_ongoing98*TBR_NSAdummy
             + numdyads +rebcap+ blood+ Home_v2svstterr+v2x_rule_Home+intensitylevel +dyad_yearatwarln 
                      +civ_deaths_dummy 
                      , data=df2, family = binomial)


# Compute robust standard errors and store the coeftest results
xf1A_98_w_controls_exceptPKO_AddrebOSV_cluster <- coeftest(
  xf1A_98_w_controls_exceptPKO_AddrebOSV, 
  vcov = sandwich::vcovCL(xf1A_98_w_controls_exceptPKO_AddrebOSV, cluster = df2$side_a_id, type = "HC0")
)
xf1A_98_w_controls_exceptPKO_AddrebOSV_cluster 

## ############AX TBR_NSAdummy  X 1998 ###################
xf1A_98_w_controls_AddrebOSV <- glm(hram ~ war_ongoing98*TBR_NSAdummy
                 + numdyads +rebcap+ blood
                 + Home_v2svstterr+v2x_rule_Home
                 +intensitylevel +dyad_yearatwarln 
                 +PKO_dummy_Home +civ_deaths_dummy 
                      , data=df2, family = binomial)

# Compute robust standard errors and store the coeftest results
xf1A_98_w_controls_AddrebOSV_cluster <- coeftest(
  xf1A_98_w_controls_AddrebOSV, 
  vcov = sandwich::vcovCL(xf1A_98_w_controls_AddrebOSV, cluster = df2$side_a_id, type = "HC0")
)
xf1A_98_w_controls_AddrebOSV_cluster 
```

```{r H1: Tables and Figures}
# Table 4: Transborder Status and Amnesty for Serious Crimes: Pre- and Post-1998
texreg(list(xf1A_98_w_controls_pre89df_AddrebOSV_cluster,  #pre98
            xf1A_98_w_controls_post89df_AddrebOSV_cluster, #post98
            xf1A_98_WO_controls_cluster,  #all w/o controls
            xf1A_98_w_controls_AddrebOSV_cluster),  #all w/ all controls
       include.ci = FALSE,
       stars = c(0.1, 0.05, 0.01, 0.001))

# For number of observation and AIC info, see this (before clustering SE): 
texreg(list(xf1A_98_w_controls_pre89df_AddrebOSV,
            xf1A_98_w_controls_post89df_AddrebOSV,
            xf1A_98_WO_controls, 
            xf1A_98_w_controls_AddrebOSV), 
       include.ci = FALSE,
       stars = c(0.1, 0.05, 0.01, 0.001))

################################################  
######## Figure 3: Pred. Probability of Amnesty by Rebel Type and Time Period  
################################################

#Predicted Probability
predDF_AX1_AddrebOSV <- ggpredict(xf1A_98_w_controls_AddrebOSV, 
                                  terms = c("war_ongoing98", "TBR_NSAdummy"),
                                  vcov.fun = "vcovCL",
                                  vcov.args = list(cluster = xf1A_98_w_controls_AddrebOSV$model$side_a_id))  


predDF_AX1_AddrebOSV$x <- factor(predDF_AX1_AddrebOSV$x, 
                                 labels = c("Pre-98", "Post-98"))
predDF_AX1_AddrebOSV$group <- factor(predDF_AX1_AddrebOSV$group, 
                                     labels = c("Local Rebel Groups", "TBRs"))

#Plotting
AX_df2 <- ggplot(predDF_AX1_AddrebOSV, aes(x = x, 
                                           y = predicted, 
                                           ymin = conf.low, 
                                           ymax = conf.high, 
                                           group = group, color = group)) +
  # dot and error bars
  geom_point(position = position_dodge(0.4), size = 3) +
  geom_errorbar(width = 0.1, position = position_dodge(0.4)) +
  labs(x = "Conflict Period",
       y = "Predicted Probability of Amnesty",
       color = "Rebel Type",
       #title = "Predicted Probabilities of Amnesty by Rebel Type and Era"
       ) +
  scale_color_manual(values = c("gray60", "black")) +
  theme_bw() + 
  theme(legend.position = "bottom")

AX_df2 #figure 2


# Exact Numbers of Pred prob for paper
print(predDF_AX1_AddrebOSV)
predDF_AX1_AddrebOSV[, c("x", "group", "predicted")]

####################################################################
### Figure 2: Marginal Effect of the 1998 Threshold on Amnesty Probability
##################################################################

# Load necessary libraries
library(marginaleffects)
library(ggplot2)

# 1. Calculate Marginal Effects (Change in probability from Pre to Post). This computes the contrast (Post-98 minus Pre-98) for each Rebel Type.

mfx_data <- avg_comparisons(
  xf1A_98_w_controls_AddrebOSV,
  variables = "war_ongoing98",
  by = "TBR_NSAdummy",
  vcov = ~side_a_id 
)

# 2. Prepare labels for visualization
mfx_data$TBR_NSAdummy <- factor(mfx_data$TBR_NSAdummy, 
                                labels = c("Local Rebel Groups", "TBRs"))

# 3. Plotting the Contrasts (Marginal Effect Plot)
mfx_plot <- ggplot(mfx_data, aes(x = TBR_NSAdummy, y = estimate, 
                                 ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") + 
  geom_point(size = 4) +
  geom_errorbar(width = 0.2) +
  labs(
    x = "Rebel Type",
    y = "Change in Predicted Probability",
  ) +
  theme_bw() +
  theme(plot.title = element_text(face = "bold"))

# Display the plot
print(mfx_plot)

# 4. Extract exact numerical results (discussed in paper)
print(mfx_data[, c("TBR_NSAdummy", "estimate", "conf.low", "conf.high", "p.value")])
```

```{r H1: Appendix A.1. Conflict After 1998}
dfbegin_post98war<- df2%>%filter(start_post1998==1)

############# Post98: AX TBR_NSAdummy  ########################
xf1A_98_w_controls_dfbegin_post98war <- glm(hram ~ TBR_NSAdummy
                   + numdyads +rebcap+ blood+ Home_v2svstterr+v2x_rule_Home+intensitylevel +dyad_yearatwarln+ +PKO_dummy_Home
                      , data=dfbegin_post98war, family = binomial)

# Compute robust standard errors and store the coeftest results
xf1A_98_w_controls_dfbegin_post98war_cluster <- coeftest(
  xf1A_98_w_controls_dfbegin_post98war, 
  vcov = sandwich::vcovCL(xf1A_98_w_controls_dfbegin_post98war, cluster = dfbegin_post98war$side_a_id, type = "HC0")
)


# Using texreg to display the models
texreg(xf1A_98_w_controls_dfbegin_post98war_cluster, 
       include.ci = FALSE,
       stars = c(0.1, 0.05, 0.01, 0.001))

# For N, AIC info, see this: 
texreg(xf1A_98_w_controls_dfbegin_post98war, 
       include.ci = FALSE,
       stars = c(0.1, 0.05, 0.01, 0.001))
```

```{r H1: Appendix A.2. The 2002 Cutoff}
################################################################
################# Split Models (M1', M2') ##########
################################################################
dfpre2002war <- df2 %>%filter(end_pre2002==1)
dfpost2002war <- df2 %>%filter(end_pre2002==0)


###### M1': Pre2002 
xf1A_pre2002 <- glm(hram ~ TBR_NSAdummy
                    + numdyads +rebcap+ blood+ Home_v2svstterr+v2x_rule_Home+intensitylevel +dyad_yearatwarln  
                    +civ_deaths_dummy
                    +PKO_dummy_Home
                        , data=dfpre2002war, family = binomial)

# Compute robust SE and store the coeftest results
xf1A_pre2002_cluster <- coeftest(
  xf1A_pre2002, 
  vcov = sandwich::vcovCL(xf1A_pre2002, cluster = dfpre2002war$side_a, type = "HC0")
)

xf1A_pre2002_cluster

###### M2': Post2002 
xf1A_post2002 <- glm(hram ~ TBR_NSAdummy
                     + numdyads +rebcap+ blood+ Home_v2svstterr+v2x_rule_Home+intensitylevel
                     +dyad_yearatwarln 
                      +civ_deaths_dummy      
                     +PKO_dummy_Home          
                     , data=dfpost2002war, family = binomial)

# Compute robust standard errors and store the coeftest results
xf1A_post2002_cluster <- coeftest(
  xf1A_post2002, 
  vcov = sandwich::vcovCL(xf1A_post2002, cluster = xf1A_post2002$side_a, type = "HC0")
)
xf1A_post2002_cluster

################################################################
################# Full Models (M3', M4') ##########
################################################################

##M3' : no control
xf1A_2002_WO_controls <- glm(hram ~ war_ongoing2002*TBR_NSAdummy, data=df2, family = binomial)

# Compute robust standard errors and store the coeftest results
Dyad_AX_WO_2002_cluster <- coeftest(
  xf1A_2002_WO_controls, 
  vcov = sandwich::vcovCL(xf1A_2002_WO_controls, cluster = df2$side_a_id, type = "HC0")
)
Dyad_AX_WO_2002_cluster

##M4': With controls  
xf1A_2002_w_controls_pko <- glm(hram ~ war_ongoing2002*TBR_NSAdummy
                    + numdyads +rebcap+ blood+ Home_v2svstterr+v2x_rule_Home+intensitylevel +dyad_yearatwarln
                     + civ_deaths_dummy
                      +PKO_dummy_Home, data=df2, family = binomial)

# Compute robust standard errors and store the coeftest results
Dyad_AX_2002_pko_cluster <- coeftest(
  xf1A_2002_w_controls_pko, 
  vcov = sandwich::vcovCL(xf1A_2002_w_controls_pko, cluster = df2$side_a_id, type = "HC0")
)
Dyad_AX_2002_pko_cluster

#####################################
############# Regression Table (Table A2: Transborder Status and Amnesty for Serious Crimes: Pre- and Post-2002)
#####################################

# Using texreg to display the models
texreg(list(xf1A_pre2002_cluster, # pre2002
  xf1A_post2002_cluster, # post2002
  Dyad_AX_WO_2002_cluster, 
  Dyad_AX_2002_pko_cluster), 
       include.ci = FALSE,
       stars = c(0.1, 0.05, 0.01, 0.001))

# For numbers, AIC see this: 
texreg(list(xf1A_pre2002,
xf1A_post2002,
xf1A_2002_WO_controls,
  xf1A_2002_w_controls_pko), 
       include.ci = FALSE,
       stars = c(0.1, 0.05, 0.01, 0.001))

#####################################
############# FIGURE A1: Predicted Probability of Amnesty by Rebel Type and Time Period (2002)
#####################################
predDF_AX1_2002 <- ggpredict(xf1A_2002_w_controls_pko, c("war_ongoing2002", "TBR_NSAdummy"),
                    vcov.fun = "vcovCL",
                    vcov.type = "HC1",
                    vcov.args = list(cluster = df2$side_a_id))


# Create condition column manually
predDF_AX1_2002$post2002war <- factor(rep(c("Pre-2002", "Post-2002"), each = 2))

# Modify the predDF to have proper labels
predDF_AX1_2002$group <- factor(predDF_AX1_2002$group, labels = c("Local Rebel Groups", "TBRs"))


library(ggplot2)
predDF_AX1_2002$post2002war <- factor(predDF_AX1_2002$post2002war, levels = c("Pre-2002", "Post-2002"))

predDF_AX1_2002$x <- factor(predDF_AX1_2002$x, labels = c("Pre-2002", "Post-2002"))
predDF_AX1_2002$group <- factor(predDF_AX1_2002$group, labels = c("Local Rebels", "TBRs"))

## plot
AX_df2_2002 <- ggplot(predDF_AX1_2002, aes(x = x, y = predicted, 
                                           ymin = conf.low, ymax = conf.high, 
                                           group = group, color = group)) +
  geom_errorbar(width = 0.1, size = 0.7, 
                position = position_dodge(0.4)) +
  geom_point(size = 3.5, position = position_dodge(0.4)) +
  
  labs(x = "Conflict Period", 
       y = "Predicted Probability of Amnesty", 
       color = "Rebel Type",
     ) +
  scale_color_manual(values = c("gray60", "black")) +
  
  theme_bw() + 
  theme(legend.position = "bottom",
        panel.grid.minor = element_blank(),
        axis.title = element_text(size = 11),
        legend.title = element_text(size = 10))

print(AX_df2_2002) 
```

```{r H1: Appendix A.3 Three Interaction}
# M1 (in A.3.)
xf1A_98_WO_controls_OSV_three <- glm(hram ~ war_ongoing98*TBR_NSAdummy*log_civ_deaths
                   , data=df2, family = binomial)

xf1A_98_WO_controls_OSV_cluster_three <- coeftest(
  xf1A_98_WO_controls_OSV_three, 
  vcov = sandwich::vcovCL(xf1A_98_WO_controls_OSV_three, cluster = df2$side_a_id, type = "HC0")
)
xf1A_98_WO_controls_OSV_cluster_three

# M2 (in A.3.)
xf1A_98_w_controls_exceptPKO_OSV_three <- glm(hram ~ war_ongoing98*TBR_NSAdummy *log_civ_deaths
                   + numdyads +rebcap+ blood+ Home_v2svstterr+v2x_rule_Home+intensitylevel +dyad_yearatwarln 
                      , data=df2, family = binomial)

xf1A_98_w_controls_exceptPKO_cluster_OSV_three <- coeftest(
  xf1A_98_w_controls_exceptPKO_OSV_three, 
  vcov = sandwich::vcovCL(xf1A_98_w_controls_exceptPKO_OSV_three, 
                          cluster = df2$side_a_id, type = "HC0")
)
xf1A_98_w_controls_exceptPKO_cluster_OSV_three 

# M3 (in A.3.)
xf1A_98_w_controls_OSV_three <- glm(hram ~ war_ongoing98*TBR_NSAdummy *log_civ_deaths 
            + numdyads +rebcap+ blood+ Home_v2svstterr+v2x_rule_Home+intensitylevel +dyad_yearatwarln 
                      +PKO_dummy_Home
                      , data=df2, family = binomial)

xf1A_98_w_controls_cluster_OSV_three <- coeftest(
  xf1A_98_w_controls_OSV_three, 
  vcov = sandwich::vcovCL(xf1A_98_w_controls_OSV_three, cluster = df2$side_a_id, type = "HC0")
)
xf1A_98_w_controls_cluster_OSV_three 

############# Regression Table ###########

# Using texreg to display the models
texreg(list(xf1A_98_WO_controls_OSV_cluster_three, 
            xf1A_98_w_controls_exceptPKO_cluster_OSV_three, 
            xf1A_98_w_controls_cluster_OSV_three), 
       include.ci = FALSE,
       stars = c(0.1, 0.05, 0.01, 0.001))

# for number of observation and AIC info, see this:  
texreg(list(xf1A_98_WO_controls_OSV_three, xf1A_98_w_controls_exceptPKO_OSV_three,xf1A_98_w_controls_OSV_three), 
       include.ci = FALSE,
       stars = c(0.1, 0.05, 0.01, 0.001))
```
